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Abstract 

We propose a practicable method for describing linear dynamics of different finite Fermi sys- 
tems. The method is based on a general self-consistent procedure for factorization of the two-body 
residual interaction. It is relevant for diverse density- and current-dependent functionals and, in 
fact, represents the self-consistent separable random-phase approximation (RPA), hence the name 
SRPA. SRPA allows to avoid diagonalization of high-rank RPA matrices and thus dwarfs the calcu- 
lation expense. Besides, SRPA expressions have a transparent analytical form and so the method 
is very convenient for the analysis and treatment of the obtained results. SRPA demonstrates high 
numerical accuracy. It is very general and can be applied to diverse systems. Two very different 
cases, the Kohn-Sham functional for atomic clusters and Skyrme functional for atomic nuclei, are 
considered in detail as particular examples. SRPA treats both time-even and time-odd dynamical 
variables and, in this connection, we discuss the origin and properties of time-odd currents and 
densities in initial functionals. Finally, SRPA is compared with other self-consistent approaches 
for the excited states, including the coupled-cluster method. 
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I. INTRODUCTION 



The time-dependent local-density-approximation theory (TDLDA) is widely used for 
description of dynamics of diverse quantum systems such as atomic nuclei, atoms and 
molecules, atomic clusters, etc. (see for more details However, even in 

the linear regime, this theory is plagued by dealing with high-rank matrices which make the 
computational effort too expensive. This is especially the case for non-spherical systems with 
their demanding configuration space. For example, in the Random Phase Approximation 
(RPA), a typical TDLDA theory for linear dynamics, the rank of the matrices is determined 
by the size of the particle-hole lph space which becomes really huge for deformed and heavy 
spherical systems. The simplest RPA versions, like the sum rule approach and local RPA 
(see hierarchy of RPA methods in p) deal with a few collective variables instead of a full 
lph space and thus avoid the problem of high-rank matrices. But these versions cannot 
properly describe gross-structure of collective modes and the related property of the Landau 
damping (dissipation of the collective motion over nearby lph excitations). 

In this connection, we propose a method llCl ] which combines accuracy and 

power of involved RPA versions with simplicity and physical transparency of the simplest 
ones and thus is a good compromise between these two extremes. The method is based 
on the self-consistent separable approximation for the two-body residual interaction which 
is factorized into a sum of weighted products of one-body operators. Hence the method is 
called as separable RPA (SRPA). It should be emphasized that the factorization is self- 
consistent and thus does not result in any additional parameters. Expressions for the one- 
body operators and their weights are inambiguously derived from the initial functional. The 
factorization has the advantage to shrink dramatically the rank of RPA matrix (usually 
from r = 10 3 — 10 6 to r = 2 — 14) and thus to minimize the calculation expense. Rank 
of SRPA matrix is determined by the number of the separable terms in the expansion for 
the two-body interaction. Due to effective self-consistent procedure, usually a few separable 
terms (or even one term) are enough for a good accuracy. Ability of SRPA to minimize the 
computational effort becomes really decisive in the case of non-spherical systems with its 
huge lph configuration space. SRPA formalism is quite simple and physically transparent, 
which makes the method very convenient for the analysis and treatment of the numerical 
results. Being self-consistent, SRPA allows to extract spurious admixtures connected with 
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violation of the translational or rotational invariance. As is shown below, SRPA exhibits 
accuracy of most involved RPA versions but for the much less expense. Since SRPA exploits 
the full lph space, it equally well treats collective and non-collective states and, what is 
very important, fully describes the Landau damping, one of the most important properties 
of collective motion. SRPA is quite general and can be applied to diverse finite Fermi 
systems (and thus to different functionals), including those tackling both time-even densities 
and time-odd currents. The latter is important not only for nuclear Skyrme functionals 
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which exploits a variety of time-even and time-odd variables but also for electronic 
functionals whose generalized versions deal with basic current densities (see e.g. jl^t). 

SRPA has been already applied for atomic nuclei and clusters, both spherical and de- 
formed. To study dynamics of valence electrons in atomic clusters, the Konh-Sham func- 
tional [hi list was exploited |7j£, 



pseudo-Hamiltonian schemes 



16 



17], in some cases together with pseudopotential and 



|. Excellent agreement with the experimental data [l8( for 
the dipole plasmon was obtained. Quite recently SRPA was used to demonstrate a non- 
trivial interplay between Landau fragmentation, deformation splitting and shape isomers in 
forming a profile of the dipole plasmon in deformed clusters |l7| . 



In atomic nuclei, SRPA was derived 19 







1911 for the demanding Skyrme functional in- 



volving a variety of densities and currents (see [20| for the recent review on Skyrme forces). 
SRPA calculations for isoscalar and isovector giant resonances (nuclear counterparts of elec- 
tronic plasmons) in doubly magic nuclei demonstrated high accuracy of the method [lol |. 

In the present paper, we give a detail, maybe even tutorial, description of SRPA, consider 
and discuss its most important particular cases and compare it with alternative approaches, 
including the equation-of-motion method in the coupled-cluster theory. We thus pursue the 
aim to advocate SRPA for researchers from other areas, e.g. from the quantum chemistry. 

The paper is organized as follows. In Section |Hj derivation of the the SRPA formalism 
is done. Relations of SRPA with other alternative approaches are commented. In Sec. 3, 
the method to calculate SRPA strength function (counterpart of the linear response theory) 
is outlined. In Section ITVl the particular SRPA versions for the electronic Kohn-Sham 
and nuclear Skyrme functionals are specified and the origin and role of time-odd currents 
in functionals are scrutinized. In Sec. 13 the practical SRPA realization is discussed. 
Some examples demonstrating accuracy of the method in atomic clusters and nuclei are 
presented. The summary is done in Sec. IVII In Appendix A, densities and currents for 
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Skyrme functional are listed. In Appendix B, the optimal ways to calculate SRPA basic 
values are discussed. 



II. BASIC SRPA EQUATIONS 

RPA problem becomes much simpler if the residual two-body interaction is factorized 
(reduced to a separable form) 

K 

^2 < h 2 P2\V res \pih 1 > a^a^ 2 a h2 a hl -> ^ [ K kk>X k X k , + r} W Y k Y k r] (1) 

h\,h,2,pi,P2 k,k'=l 

where 

X k = ^< p\X k \h > a+a h , Y k = ^< p\Y k \h > a+a h 

ph ph 

are time-even and time-odd one-body operators, respectively. Further, a+ (a^) is the creation 
(annihilation) operator of the particle state p (hole state h); K is the number of the separable 
terms. 

Conceptually, the self-consistent procedure outlined below was first proposed in j^lj ]. 
A. Time-dependent Hamiltonian 

The system is assumed to undergo small-amplitude harmonic vibrations around HF 
ground state. The starting point is a general time-dependent energy functional 

E(J a (f,t)) = J H(J a (r,t))dr (2) 

depending on an arbitrary set of densities and currents defined through the corresponding 
operators as 

occ 

J a (f,t) =< ¥(f)|J a (fO|*(f) >=J2^h(r,t)J a (^ h (r,t) (3) 

h 

where \l/(t) is the many-body function of the system as a Slater determinant, and tp* h is the 
wave function of the hole (occupied) single-particle state. In general, the set (jHJ) includes 
both time-even and time-odd densities and currents, see examples in the Appendix A. 
Time-dependent mean- field Hamiltonian directly follows from 0-©: 



In the small-amplitude regime, the densities are decomposed into static part and small 
time-dependent variation 

J a (f,t) = J a (r) + 5J a (f,t). (5) 

Then, to the linear order for 5J a (r,t), the mean-field Hamiltonian (@J) can be decomposed 
into static and time-dependent response parts 

h(r,t) = h {r) + h res (f,t), 

= J2^=jUr) + J2^s^h=jSJAr,t)J a (^ (6) 

a ' a,a' 

and thus we get the time-dependent Hamiltonian h res (f,t) responsible for the collective 
motion. 

B. Scaling perturbation 

Now we should specify the response Hamiltonian h res (f,t). For this aim, we use the 
scaling transformation and define the perturbed many-body wave function of the system as 

K 

Mt)>= Y[exp[-iq k (t)P k ]exp[-i Pk (t)Q k }\0 > . (7) 
fe=i 

Here both the perturbed wave function \^f(t) > and static ground state wave function |0 > are 
Slater determinants; Q k {r) and P k {r) are generalized coordinate (time-even) and momentum 
(time-odd) hermitian operators with the propoerties. 

Qk = Qti TQ k T 1 = Q k , 
P k = i[H, Q k ] ph = P+, TPkf- 1 = -P k , (8) 

They generate T-even and T-odd harmonic deformations q k (t) and p k (t); T is the time 
inversion operator. 

Using Eqs. Q and (J7J), the transition densities read 

5J a (f,t) =< tf(t)|J a |tf(t) > - < 0|J a |0 >= (9) 
= *Z){9*(*) < 0\[P k , J a (r)]\0 > +Pk(t) < 0|[Qjfc, J a (f)]|0 >} 

k 

and the response Hamiltonian (JHJ) is 

h(r,t) = Y,Mt)Mr)+Pk(t)Y k (r)} (10) 

sk 
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where all the r-dependent terms are collected into the hermitian one-body operators 

= l Y}jJJ^)j=3 < °l[A,i Q ']|0 > J Q (r), (11) 

act' 



aa' 



with the properties 

■X-k = ■> TX k T 1 = X k , X k = X k , (13) 
Y k = Y k + , TY k T- l = -Y k , Y k * = -Y k . (14) 

As is shown below, X k and Y k are just the T-even and T-odd operators to be exploited in 
the separable expansion (JTJ). 

In the derivation above, we used the property of hermitian operators with a definite 
T-parity 

< 0|[i,£]|0>= 0, if TAT' 1 = TBT~ 1 (15) 

which states that the average of the commutator vanishes if the operators A and B are of 
the same T-parity. This property allows to classify operators with a definite T-parity in the 
SRPA formalism and thus to make the formalism simple and transparent. For example, in 
Eqs. (fTT j) -(fT2 j) . T-even and T-odd densities J a (r) contribute separately to X k and Y k . 

To complete the construction of the separable expansion $IJ , we should yet determine the 
strength matrices K kk < and i] kk r. This can be done through variations of the basic operators 

5X k (t) =< tt(t)|X fc |tf(t) > - < 0|X fe |0 >= (16) 
= z < 0\[P k ',X k }\0 >=-Y^qk'{t)K- k } k , 

k' k' 

5Y k (t) =< *(t)|y fc |*(t) > - < 0\Y k \0 >= (17) 



= tj2pk'(t) < o|[4',n]|o >= - 

y k' 
where we introduce symmetric inverse strength matrices 

«vl = «S = -*< 0|[A', |0 >= (18) 
= / drY}^^}<^\\hJa\\Q><^\\Pk'JA^>^ 



aa 
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%! = ^ = -^<o|[Q fc ',n]|o> (19) 

= / dfY j [-^-}<0\[Q k ,J a ]\0><0\[Q kl ,J a ,]\0> . 
J , 0J a/ dJ a 

Then one gets 

-J2 K k>kSx k (t) = q k >(t), (20) 

k 

-J2^'kSY k (t) = Pk'(t) (21) 

sk 

and the response Hamiltonian ()10|) acquires the form 

fc(r,t) = - J2{ K kk'SX k (t)Xh>(r) + r/ fcfc ^f fe (t)F F (f)}. (22) 

Following |l|, the response Hamiltonian (}2*2*|) leads to the same eigenvalue problem as the 
separable Hamiltonian 

Hrpa = h + V res , (23) 

with 

Vres = — - ^^[^fcfc'^fc^fc' + llkk'Y k Y k i] (24) 

(see also discussion in the next subsections). 

In principle, we already have in our disposal the SRPA formalism for description of the 
collective motion in space of collective variables. Indeed, Eqs. ffTTj) . (JT2J), (fTHj). and (fT9^) 
deliver one-body operators and strength matrices we need for the separable expansion of the 
two-body interaction. The number K of the collective variables q k {t) and p k (t) and separable 
terms depends on how precisely we want to describe the collecive motion (see discussion in 
Section 4). For K — 1, SRPA converges to the sum rule approach with a one collective mode 

n 

For K > 1, we have a system of K coupled oscillators and SRPA is reduced to the local 

nn 

RPA [6, 24] suitable for a rough description of several modes and or main gross-structure 
efects. However, SRPA is still not ready to describe the Landau fragmentation. For this 
aim, we should consider the detailed lph space. This will be done in the next subsection. 



C. Introduction of lph space 



Collective modes can be viewed as superpositions of lph configurations. It is convenient to 
define this relation by using the Thouless theorem which establishes the connection between 
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two arbitrary Slater determinants [25]. Then, the perturbed many-body wave function reads 

|^(t)>=(l + ^ Cp ,(t)i+)|* > (25) 

ph 

where 

A; h = a\a h (26) 

is the creation operator of lph pair and 

c P h(t) = c + ph e*« + ctfT™ (27) 

is the harmonic time-dependent particle- hole amplitude. T-even quit) and T-odd Pk(t) col- 
lective variables can be also specified as harmonic oscillations 

q k {t) = q k cos{uot) = ^q k (e iujt + e~ iujt ), 

p k {t) = p k sin{uit) = ^Pkie™ - e"*"*). (28) 

Substituting (JTUJl and (|2*5j) into the time-dependent HF equation 

ij t \^{t) >= (k + hre.(t))\V(t) >, (29) 

one gets, in the linear approximation, the relation between c ph and collective deformations 
q k and p k 

± 1 Efc'fe' < ph\X k >\0 > Tipk' < ph\Y k >\0 >] 

Cph ~ 2 T^o ' (30) 

where e p h is the energy of lph pair. 

In addition to Eqs. (fT7 |) -([T% |) . the variations SX k (t) and 5Y k (t) can be now obtained with 
the alternative perturbed wave function (|25j) : 

5X k ,{t)=Y,{c P h{ty <ph\X k ,\0> +c ph (t) <0\X k ,\ph>), (31) 

ph 

6Y k ,(t)=Y,{c ph {ty <ph\Y k ,\0> +c s ph (t) <0\Y k ,\ph>). (32) 

ph 

It is natural to equate the dynamical variations of the basic operators 5X k and 5Y k , ob- 
tained with the scaling (JJJ) and Thouless f)25j) perturbed wave functions. This provides the 
additional relation between the amplitudes c^ h and deformations q k and p k and finally result 
in the system of equations for the unknowns q k and p k . 
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By equating ([T7|l - (P|l and we get 

- ^ ?*(*)«**' = E^W* <P /i l^'l°> + c ^(*) <0|X fc ,|p/i>), (33) 

k ph 

-^^(t)^" 1 , = ^(^(t)* <ph\Y k ,\0> +c ph (t) <0\Y k ,\ph>). (34) 

A; ph 

Substituting (|?7j) - (|Snj) into these expressions and collecting, for example, the terms at e ZUJt , 
one finally gets 

E - ^}+PkF ( k X k Y) } = 0, 



with 



4^ =E 2 Z 2 {<P^|0>*<pfe|^|0>(gph + a;) (36) 



ph £ ph~ u 



+ < p/i|X fc |0 >< 0\X k >\ph > (e ph - u)}, 
=~»E 2 _ 2 {<ph\Y k \0>*<ph\X k ,\0>(e ph + uj) (37) 



+ < p/i|Y fc |0 >< 0\X k ,\ph > {e ph -uj)}, 
F^k Y) =^E 2 Z 2 {<P^|0>*<pfe|n'|0>(g P / 1 + ^) (38) 



+ < ph\X k \h >< 0\Y k ,\ph > (e ph - u)}, 
F kl Y) =E 2 Z 2 {<pfe|n|0>*<P^fe|0>(g P ft + ^) (39) 



+ < p/i|r fc |o >< o|y fc /|p/i > - w)}. 



Equating determinant of the system to zero, we get the dispersion equation for RPA 
eigenvalues uj v . 

D. Normalization condition 

By definition, RPA operators of excited one-phonon states read 

Qt = \ - C h A ph ] (40) 

ph 



and fulfill 

[Qu,Qt>] = 8*y, lQt,Qt>\ = [Qu,QA = o, (41) 

where A^ h and are given by (J2fij) and (jSOjl , respectively. In the quasiboson approximation 
for Ap h , the normalization condition [QuiQt] = 1 results in the relation 

£{(<£) a - (#) 2 } = 2 - (42) 

Using (|30p. it can be reformulated in terms of the RPA matrix coefficients (J5fi)l - (J59Jl : 

£{(<&) 2 - (St) 2 ) (43) 

du)v +^k>P k du)v +Pk>P k du)v \-lN v . 

kk' 

The variables q k and p v k should be finally normalized by the factor 1 / y/N v . 



E. General discussion 



Eqs. (dH, dH, (ffH), dH, d30|), and g^-© constitute the basic SRPA for- 

malism. It is worth now to comment some essential points: 

• One may show (e.g. by using a standard derivation of the matrix RPA) that the separa- 
ble Hamiltonian (J^H-flHD with (gUJ) results in the SRPA equations (jH^-flU if to express 
unknowns c£j~ through q k and p k . Generally, RPA equations for unknowns c£j~ require the 
RPA matrix of the high rank equal to size of the lph basis. The separable approximation 
allows to reformulate the RPA problem in terms of much more compact unknows q% and 
p k (see relation (jHO)) and thus to minimize the computational effort. As is seen from (|35jl . 
the rank of the SRPA matrix is equal to a double number K of the separable operators and 
hence is low. 

• The number of RPA eigen-states u is equal to the number of the relevant lph configura- 
tions used in the calculations. In heavy nuclei and atomic clusters, this number ranges the 
interval 10 3 -10 6 . For every RPA state v, Eq. (}35|) delivers a particular set of the amplitudes 
q"k and p u sk which, following Eq. (|HUj) . self-consistently regulate relative contrubutions of 
different T-even and T-odd oscillating densities to the i/- state. 

• Eqs. dHJ), (H2); JIBI)> (EH) relate the basic SRPA values with the starting functional and 
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input operators Q k and Pk by a simple and physically transparent way. This makes SRPA 

very convenient for the analysis and treatment of the obtained results. 

• It is instructive to express the basic SRPA operators via the separable residual interaction 



Xk — [V res , Pk]ph, Y k — [V res , Qk\ph (44) 

where the index ph means the lph part of the operator. It is seen that the T-odd operator 
Pk retains the T-even part of V res to build Xf~. Vice versa, the commutator with the T-even 
operator Qk keeps the T-odd part of V res to build Yj,. 

• Some of the SRPA values read as averaged commutators between T-odd and T-even 
operators. This allows to establish useful relations with other models. For example, (|18jl. 
(HH and flSJ) give 

%l = ~i< 0\[P k/ ,X k }\0 >= -i < Q\[P k ,, [V res ,P k }}\0 >, (45) 

Vk'l = ~ % < 0|[4',^]|0 >= -Z < 0|[4', [Vres,Qk]]\0 > . (46) 

The similar double commutators but with the full Hamiltonian (instead of the residual 
interaction) correspond to m 3 and mj sum rules, respectively, and so represent the spring 
and inertia parameters [24J in the basis of collective generators Qk and A. This allows to 



establish the connection of the SRPA with the sum rule approach 



24 1 . 



22, 231 and local RPA 



Besides, the commutator form of SRPA values can considerably simplify their calculation 
(see discussion in the Appendix B). 

• SRPA restores the conservation laws (e.g. translational invariance) violated by the static 
mean field. Indeed, let's assume a symmetry mode with the generator P sym . Then, to keep 
the conservation law [H, -P sym ] = 0, we simply have to include P sym into the set of the input 
generators P k together with its complement Q sym = i[H,P sym \. 

• SRPA equations are very general and can be applied to diverse systems (atomic nuclei, 
atomic clusters, etc.) described by density and current- dependent functionals. Even Bose 
systems can be covered if to redefine the many-body wave function (J25)) exhibiting the per- 
turbation through the elementary excitations. In this case, the Slater deterninant for lph 
excitations should be replaced by a perturbed many-body function in terms of elementary 
bosonic excitations. 

• In fact, SRPA is the first TDLDA iteration with the initial wave function (|7jl. A single 
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iteration is generally not enough to get the complete convergence of TDLDA results. How- 
ever, SRPA calculations demonstrate that high accuracy can be achieved even in this case if 
to ensure the optimal choice of the input operators Qk and Pk and keep sufficient amount of 
the separable terms (see discussion in Sec. 5). In this case, the first iteration already gives 
quite accurate results. 



• There are some alternative RPA schemes also de 
the two-body residual interaction, see e.g. [21 



ivering self-consistent factorization o 



28, 



291 ] for atomic nuclei and 
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for atomic clusters. However, these schemes are usually not sufficiently general. Some of 
them are limited to analytic or simple numerical estimates next ones start from 

phenomenological single-particle potentials and thus are not fully self-consistent [28|, the 
others need a large number of the separable terms to get an appropriate numerical accuracy 
(29, 31]. SRPA has evident advantages as compared with these schemes. 
• After solution of the SRPA problem, the Hamiltonian ()23j) - ([2*4*j) is reduced to a composition 



of one-phonon RPA excitations 



(47) 



where one-phonon operators are given by (j4U|) — ()41|) . Then, it is easy to get expressions of 
the equation-of-motion (EOM) method: 



(48) 



So, SRPA and EOM with the Hamiltonian (|23j ) -([24 p are equivalent. This allows to to 
establish the connection between the SRPA and couled-cluster EOM method with the single 



reference (see for reviews |3z 
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34j, 3)- SRPA 

uses the excitation operators involving 
only singles (lph) and so generally carries less correlations than EOM-CC. At the same 
time, SRPA delivers very elegant and physically transparent calculations scheme and, as is 
shown in our calculations, the correlations included to the SRPA are often quite enough 
to describe linear dynamics. It would be interesting to construct the approach combining 
advantages of SRPA and EOM-CC. 



III. STRENGTH FUNCTION 



In study of response of a system to external fields, we are usually interested in the 
average strength function instead of the responses of particular RPA states. For example, 
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giant resonances in heavy nuclei are formed by thousands of RPA states whose contributions 
in any case cannot be distinguished experimentally. In this case, it is reasonable to consider 
the averaged response described by the strength function. Besides, the calculation of the 
strength function is usually much easier. 

For electric external fields of multipolarity EXfi, the strength function can be defined as 

Sl{E\(j,;u) = 5>, L M A %C(^-^) (49) 

where 

C( "-^ ) = i(c-^A(A/2) 2 (50) 
is Lorentz weight with an averaging parameter A and 

M ^ = \Y.< p h \f^\° > K' h + c vi) ( 51 ) 

ph 

is the transition matrix element for the external field 

Am = TT W(Y AM + y A tj. (52) 

It is worth noting that, unlike the standard definition of the strength function with using 
8 (id — uj u ), we exploit here the Lorentz weight. It is very convenient to simulate various 
smoothing effects. 

The explicite expression for (J49j) can be obtained by using the Cauchy residue theorem. 
For this aim, the strength function is recasted as a sum of v residues for the poles z = 
±c<v Since the sum of all the residues (covering all the poles) is zero, the residues with 
z = ±u u (whose calculation is time consuming) can be replaced by the sum of residies with 
z = 10 ± i(A/2) and z = ±e p h whose calculation is much less expensive (see details of the 
derivation in [q]). 

Finally, the strength function for L=0 and 1 reads as 



Sl(E\[i, uo) = — 3 

7T 



z L det \B(z) 
detlFMI 



(53) 

-i(A/2) 



K p ,K h >0 

+ 2V2 J2 4 h <ph\f Xll \0> 2 ((uj-e ph ). 

ph 

The first term in (p)3*j) is contribution of the residual two-body interaction while the second 
term is the unperturbed (purely lph) strength function. Further, F(z) is determinant of the 
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RPA symmetric matrix (|35|) of the rank 2K, where K is the number of the initial operators 
Qk- The symmetric matrix B of the rank (2K + 1) is defined as 



B n n'{z) = F nn r(z), 



(54) 



B 2K +l,2K+l(z) = 0, B 2K +l,n(z) = B nt 2K+l(z) = A n (z) 



where n, n' = 1, ...,2K and left (right) indexes define lines (columns). 
The values A n (z) read 



K p ,K h >0 

A^ 2k _ 1 {EX^z) = A £ 




(55) 




They form the right and low borders of the determinant B thus fringering the RPA de- 
terminant F. The values —A n (z) in the most right column have the same indices as the 
corresponding strings of the RPA determinant. The values A n (z) in the lowest line have the 
same indices as the corresponding columns of the RPA determinant. 

Derivation of the strentth function, given above, deviates from the standard one in the 
lineary response theory. Besides, the SRPA deals with the Lorentz weight instead of 6(u—u v ) 
used in the linear response theory. At the same time, SRPA strength function and lineary 
response theory are conceptually the same approaches. Since the linear response theory is 
widely used in the coupled-cluster (CC) method, it would be interesting to consider the 
implementation of SRPA stength function method to CC. The linear response theory is 
widely used in the coupled-cluster (CC) method. In this connection, it would be interesting 
to enlarge the SRPA stength function method to CC. 

IV. PARTICULAR CASES FOR CLUSTERS AND NUCLEI 
A. Kohn-Sham functional for atomic clusters 

Kohn-Sham functional for atomic clusters reads 



E tot (t) = E kin (t) + E xc (t) + E c {t) = / drH(p(f,t)) 



(56) 
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where 



E kin (t) = ^— / dfr(f,t) , (57) 



2m e 



E xc {t) = J dfp{f, t) e xc {p (f, t)) , (58) 

Ec (t) = e -[ [ dfa ^ ^^)-^)^*)-^i)) (59) 

are kinetic, exchange-correlation (in the local density approximation), and Coulomb terms, 
respectively. Further, Pi(r) is the ionic density and p(f,t) and r(f,t) are density and kinetic 
energy density of valence electrons. 

In atomic clusters, oscillations of valence electrons are generated by time-dependent varia- 
tions of the electronic T-even density p(f, t) only. So, one may neglect in the SRPA formalism 
all T-odd densities and their variations Pk(t)- This makes SRPA equations especially simple. 
In particular, the density variation (JJJJ is reduced to 

6p(f,t) = <0| [A, p(r)]|0> (60) 

k 

K p ,K h >0 

= -4i^2q k (t) <ph\P k \0 > $l<ph\p\0 > 



ph 



= % 

k 

where 

<ph\P k \0 >= 2is ph <ph\Q k \ >, (61) 

<W) = "^(VP(r) • \/Q k {r) + 2p(r)AQ k (r)). (62) 

Here, 3?< ph\p\0 > is the transition density and p(r) is the static ground state density of 
valence electrons. 

It is seen from ()60j ) -([6i p that there are two alternative ways to calculate the density 
variation: i) through the transition density and matrix elements of Qfc-operator and ii) 
through the ground state density. The second way is the most simple. It becomes possible 
because, in atomic clusters, V res has no T-odd Y^-operators and thus the commutator of Q k 
with the full Hamiltonian is reduced to the commutator with the kinetic energy term only: 

P k = i[H, Qk] P h = i[ho, Q k ] p h = _ ^^-[V > Qk] P h- (63) 
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This drastically simplifies SRPA expressions and allows to present them in terms of the 
static ground state density. The scaling tranformation (JJJ) loses exponents with pk{t) and 
reads 

K 

!*(*)>= l[exp{-iq k (t)P k }\0 > . (64) 



k=l 

Other SRPA equations are reduced to 



5p5p 

d 2 H xc . 2 f Spin 



M?) = A—\p=P < 0|[A,p]|0 > p{r 



^^--Mrl + e'Jdr^, (65) 

*vl =K>w = -i<0\[Pk>,Xk]\0> 

= J df{^} p= ?<0\[P k ,p]\0><0\[P k , : p]\0> 

= - J drX k (r)5p(r), (66) 

E^4^-*S1} = 0. (67) 

it 

,± _ _ 1 J2 k r k <ph\X k \0> 

ph ~ 2 e ph ±u u ' m 

E< W) a - W) a > = E = 2iV - (69) 

fefc' 

It is seen that the basic operator ()65j) and strength matrix ()66)1 have now simple expressions 
via Sp(f) from (j62j) . The operator (|63|) has exchange- correlation and Coulomb terms. For 
electric multipole oscillations (dipole plasmon, ...), the Coulomb term dominates. 

In recent years, there appear some new functionals where the current of electrons in- 



stead of their density is used basic variable 



3- 



SRPA equations for this case can be 



straightforwardly obtained from the general formalism given in Sec. 2. 



B. Skyrme functional for atomic nuclei 

Nuclear interaction is very complicated and its explicit form is still unknown. So, in 



11 



practice different approximations to nuclear interaction are used. Skyrme forces 
represent one of the most successful approximations where the interaction is maximally 
simplified and, at the same time, allows to get accurate and universal description of both 
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ground state properties and dynamics of atomic nuclei (see [20( a for recent review). Skyrme 
forces are contact, i.e. ~ 5{f\ — F2), which minimizes the computational effort. In spite of 
this dramatic simplification, Skyrme forcese well reproduce properties of most spherical and 
deformed nuclei as well as characteristics of nuclear matter and neutron stars. Additional 
advantage of the Skyrme interaction is that its parameters are directly related to the basic 
nnclear properties: m _ si bi«arradii, masS e S and binding ene rgl e S , etc. SRPA 



for Skyrme forces was derived in |9j, 



10 
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26). 



Although Skyrme forces are relatively simple, they are still much more demanding than 
the Coulomb interaction. In particular, they deal with a variety of diverse densities and 



currents. The Skyrme functional reads |12 



36 



E = J dr(n kin + H S k(Ps, r s , a s ,] s , J s ) + Wc(p P )) , 



where 



Ti.sk 



2m 



-T. 



drp p (f) 



3 2-3,1 



2 / rFy ' r-fT PV 4 V 



^0 2 _ ^0 



2 h 2 



U 



(70) 

(71) 
(72) 
(73) 



+b 1 (pr-?)-b' 1 J2(psr s -J 2 s ) 

s 

-64 (pV3 + ^(Vxj))-6^ (ps(V3 s ) + a s ■ (V x I 

s 

+64 (of - 3 2 ) + \ { S 'T> 



C>2 



are kinetic, Coulomb and Skyrme terms respectively. The isospin index s = n,p covers 
neutrons (n) and protons (p). Densities without this index involve both neutrons and 
protons, e.g. p = p p + p n . Parameters b and a are fitted to describe ground state properties 
of atomic nuclei. 

The functional includes diverse densities and currents, both neutrons and protons. They 
are naturally separated into two groups: 1) T-even density p s (r), kinetic energy density 
T s (f) and spin orbital density $s s (r) and 2) T-odd spin density u s (r), current j a (r) and 
vector kinetic energy density T s (r). Explicit expressions for these densities and currents, as 
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well as for their operators, are given in the appendix A. Only T-even densities contribute to 
the ground state properties of nuclei with even numbers of protons and neutrons (and thus 
with T-even wave function of the ground state). Instead, both T-even and T-odd densities 
participate in generation of nuclear oscillations. Sec. 2 delivers the SRPA formalism for this 
general case. 



C. T-odd densities and currents 



Comparison of the Kohn-Sham and Skyrme functionals leads to a natural question why 
these two functionals exploit, for the time-dependent problem, so different sets of basic 
densities and currents? If the Kohn-Sham functional is content with one density, the Skyrme 
forces operate with a diverse set of densities and currents, both T-even and T-odd. Then, 
should we consider T-odd densities as genuine for the description of dynamics of finite many- 
body systems or they are a pequliarity of nuclear forces? This question is very nontrivial 
and still poorly studied. We present below some comments which, at least partly, clarify 
this point. 

Actual nuclear forces are of a finite range. These are, for example, Gogny forces 
representing more realistic approximation of actual nuclear forces than Skyrme approxima- 
tion. Gogny interaction has no any velocity dependence and fulfills the Galilean invariance. 
Instead, two-body Skyrme interaction depends on relative velocities k = l/2i ■ (Vi — V2), 
which just simulates the finite range effects J^. 

The static Hartree-Fock problem assumes T-reversal invariance and T-even single-particle 
density matrix. In this case, Skyrme forces can be limited by only T-even densities: p s {r), 
T s (f) and ^> s (r). In the case of dynamics, the density matrix is not already T-even and 
aquires T-odd components [12]. This fact, together with velocity dependence of the Skyrme 
interaction, results in appearance in the Skyrme functional of T-odd densities and currents: 
s s (r), j s (r) and T s (f) |l2u36|. Hence the origin of T-odd densities in the Skyrme functional. 
However, this is not the general case for nuclear forces. 

As compared with the Kohn-Sham functional for electronic systems, the nuclear Skyrme 
functional is less genuine. The main (Coulomb) interaction in the Kohn-Sham problem is 
well known and only exchange and corellations should be modeled. Instead, in the nuclear 
case, even the basic interaction is unknown and should be approximated, e.g. by the simple 
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contact interaction in Skyrme forces. 

The crudeness of Skyrme forces has certain consequences. For example, the Skyrme 
functional has no any exchange-correlation term since the relevant effects are supposed to be 
already included into numerous Skyrme fitting parameters. Besides, the Skyrme functional 
may accept a diverse set of T-even and T-dd densities and currents. One may say that 

its specific construction. 



4l|) do not exploit T-odd 



T-odd densities appear in the Skyrme functional partly because o: 
Indeed, other effective nuclear forces (Gogny j3] , Landau-Migdal 
densities and currents for description of nuclear dynamics. 

Implementation of a variety of densities and currents in the Skyrme fuctional has, how- 
ever, some advantages. It is known that different projectiles and external fields used to 
generate collective modes in nuclear reactions are often selective to particular densities and 
currents. For example, some elastic magnetic collective modes (scissors, twist) are asso- 
ciated with variations in the momentum space while keeping the common density p s (f) 
about constant. T-odd densities and currents can play here a significant role while Skyrme 
forces obtain the advantage to describe the collective motion by a natural and physically 
transparent way. 

Relative contributions of T-odd densities to a given mode should obviously depend on 
the character of this mode. Electric multipole excitations (plasmons in atomic clusters, EX 
giant resonances in atomic nuclei) are mainly provided by T-even densities (see e.g. |l9j). 
Instead, T-odd densities and currents might be important for magnetic modes and maybe 
some exotic (toroidal, ...) electric modes. 

It worth noting that T-odd denstities appear in the Skyrme functional in the specific 
combinations p s r s — jg, p s (VS s ) + a s • (V x j s ), and a s T s — 3^. Following [37(, this ensures 
Skyrme forces to fulfill the local gauge invariance (and Galilean invariance as the particular 
case). The velocity- independent finite-range Gogny forces also keep this invariance. Being 
combined into specific combinations, T-odd densities do not require any new Skyrme param- 
eters j^. So, the parameters fitted to the static nuclear properties with T-even densities 
only, are enough for description of the dynamics as well. 

38| | for electronic systems is usually imple- 



The time-dependent density functional theory 
mented at adiabatic local density approximation (ALDA) when density and single-particle 
potential are supposed to vary slowly both in time and space. Last years, the current- 
dependent Kohn-Sham functionals with a current density as a basic variable were introduced 



19 



to treat the collective motion beyond ALDA (see e.g. [13]). These functional are robust for 
a time-dependent linear response problem where the ordinary density functionals become 
strongly nonlocal. The theory is reformulated in terms of a vector potential for exchange 
and correlations, depending on the induced current density. So, T-odd variables appear in 
electronic functionals as well. 

In general, the role of T-odd variables in dynamics of finite many-body systems is still 
rather vague. This fundamental problem devotes deep and comprehensive study. 



V. CHOICE OF INITIAL OPERATORS 

It is easy to see that, after choosing the initial operators Qk{r), all other SRPA values 
can be straightforwardly determined following the steps 

Qk =*> < \[Qk,J a ]\ > => Y k , rill Pk => < \[Pk,J a ]\ > => X k , K kk \. 

As was mentioned above, the proper choice of initial operators Qfc(r) is crucial to achieve 
good convergence of the separable expansion (0) with a minimal number of separable oper- 
ators. 

SRPA itself does not give a recipe to determine Qkif)- But choice of these operators 
can be inspired by physical and computational arguments. The operators should be simple 
and universal in the sense that they can be applied equally well to all modes and excitation 
channels. The main idea is that the initial operators should result in exploration of different 
spatial regions of the system, the surface and interior. This suggests that the leading scaling 
operator should have the form of the applied external field in the long-wave approximation, 
for example, 

Qt 1 (f)=r\Y x ,(Q)+h.c). (74) 

Such a choice results in the separable operators (JTTj). (fT2|) and (pj5j) most sensitive to the 
surface of the system. This is evident in where 5p{r) oc XJp{r) is peaked at the surface. 
Many collective oscillations manifest themselves as predominantly surface modes. As a 
result, already one separable term generating by (J73)l usually delivers a quite good description 
of collective excitations like plasmons in atomic clusters and giant resonances in atomic 
nuclei. The detailed distributions depends on a subtle interplay of surface and volume 
vibrations. This can be resolved by taking into account the nuclear interior. For this aim, 
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the radial parts with lamer nnwers r^ +Tlk Y\.. and snherioal Ressel functions ran be used, much 
similar as in t 
(HQ, (H2J) and □ 
found a most 



2 -n 



Na^ 



5 2 = 0.355 




5 2 = 


0.33 


5 4 = 


0.08 




FIG. 1: Photoabsorption cross section for the dipole plasmon in axially deformed sodium clusters, 
normalized to the number of valence electrons N e . The parameters of quadrupole and hexadecapole 
deformations are given in boxes. The experimental data 39] (triangles) are compared with SRPA 



results given as bars for RPA states and as the strength function ()49|) smoothed by the Lorentz 
weight with A = 0.25 eV. Contributions to the strength function from \x =0 and 1 dipole modes 
(the latter has twice larger strength) are exhibited by dashed curves. The bars are given in eVA 2 . 



For description of the dipole plasmon in atomic clusters, the set of hermitian operators 



Qi"W = r x ^(Y Xkll (Q) + Y^JQ)) (75) 

lr^ |. As is seen from Fig. 1, we 



with XkUk = 10, 12, 14 and fj, — 0, 1 is usually enough 
successfully reproduce gross structure of the dipole plasmon in light axially-deformed sodium 
clusters (some discrepancies for the lightest cluster Na^ arise because of the roughness of 
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(0 [MeV] 



FIG. 2: Isoscalar E2 and isovector El giant resonances in 40 Ca calculated with SkM* forces. The 
results are exhibited for full (exact) RPA (solid curve) and SRPA with k = 1 (dotted curve). 

the ionic jellium approximation for smallest samples). Already one initial operator is usually 
enough to reproduce the energy of the dipole plasmon and its branches but in this case the 
plasmon acquires some artificial strength in its right flank and thus the overestimated width 
|7j. This problem can be solved by adding two more initial operators. The calculations for 
a variety of spherical alkali metal clusters [16] as well as for deformed clusters of a medium 
size show that SRPA correcly describes not only gross structure of the dipole plasmon 
but also its Landau damping and width. 

For the description of giant resonances in atomic nuclei, we used the set of initial operators 

Q 

Q k (r) = R k (r)(Y Xfl (n)+h.c.) (76) 

with 

{r\ k = l 

* M (77) 

Q\ — a k~pr~) a 2 = 0.6 ,a 3 = 0.9 ,a 4 = 1.2 

-KdifT 

where R^m is the diffraction radius of the actual nucleus and z\ is the first root in of the 
spherical Bessel function j\(z x ) = 0. The separable term with k — 1 is mainly localized at 
the nuclear surface while the next terms are localized more and more in the interior. This 
simple set seems to be a best compromise for the description of nuclear giant resonances in 
light and heavy nuclei. Fig. 2 demonstrates that already one separable term (k=l) can be 
enough to get a reasonable agreement with the exact results. For k=l, the calculations are 
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especially simple and results are easily analyzed. 

The sets (|75 ]) -(|77 j) are optimal for description of electric collective modes (EX plasmons 
in clusters and giant resonances in nuclei). For description of magnetic modes, the initial 
operator should resemble the T-odd magnetic external field. So, in this case we should start 
from the initial operators Pk in the form of the magnetic multipole transition operator in the 
long-wave approximation. The T-even operators Qk are then obtained from the connection 
Q k =i[H,P k ]. 

VI. SUMMARY 

We presented fully self-consistent separable random-phase-approximation (SRPA) 
method for description of linear dynamics of different finite Fermi-systems. The method 
is very general, physically transparent, convenient for the analysis and treatment of the 
results. SRPA drastically simplifies the calculations. It allows to get a high numerical ac- 
curacy with a minimal computational effort. The method is especially effective for systems 
with a number of particles 10 — 10 3 , where quantum-shell effects in the spectra and responses 
are significant. In such systems, the familiar macroscopic methods are too rough while the 
full-scale microscopic methods are too expensive. SRPA seems to be here the best com- 
promise between quality of the results and the computational effort. As the most involved 
methods, SRPA describes the Landau damping, one of the most important characteristics of 
the collective motion. SRPA results can be obtained in terms of both separate RPA states 
and the strength function (linear response to external fields). 

The particular SRPA versions for electronic Kohn-Sham and nuclear Skyrme functional 
were considered and examples of the calculations for the dipole plasmon in atomic clusters 
and giant resonances in atomic nuclei were presented. SRPA was compared with alternative 
methods, in particular with EOM-CC. It would be interesting to combine advantages of 
SRPA and couled-cluster approach in one powerful method. 
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Appendix A: Densities and currents for Skyrme functional 

In Skyrme forces, the complete set of the densities involves the ordinary density, kinetic- 
energy density, spin-orbital density, current density, spin density and vector kinetic-energy 
density: 

occ 

fpf- 1 = p 



frf" 1 



T^T' 1 = 3 



p 8 (f,t) = '%2<P* h (r,t)<p h (r,t), 

hes 
occ 

T S (?,t) = Y,v<p* h (f f ,tyv<p h (f,t), 

hes 

occ 

%(f,t) = -tJ2vl(r,t)(V x <?K(f,t), 

hes 
. occ 

Ur, t) = -- [p* h (f, t) V^(r, t) - V^(f, tK(f, tj\ , fjf - 1 = -j 

hes 

occ 

B s{r) =^<p* h (r,t)a<p h (r,t), 

hes 
occ 

f s (r) = J2^l(r,t)ff-V Vh (f,t), 



fat' 1 = -a 
faf- 1 = -a 



hes 



where the sum runs over the occupied (hole) single-particle states h. The associated opera- 
tors are 



Ps{r) = J^5(fi -f), 
i=\ 

i=\ 

= $^(rl -f)Vx<?, 

Js(rl = \j2{v,5(r t -r)}, 

i=i 

i=i 

Ns 



i=l 



where a is the Pauli matrix, N s is number of protons or neutrons in the nucleus. 
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Appendix B: Presentation of responses and strength matrices through the matrix 
elements 



Responses < 0\[Pk, <7<y]|0 > and < 0|[Qfc, J a ']\Q > m (fTT|) - (fT^|) and inverse strength ma- 
trices in ()18|) -()19| 1 read as the averaged commutators 

<0|[i,£]|0> with TAT' 1 = -TBT- 1 . (78) 

Calculation of these values can be greatly simplified if to express them through the lph 
matrix elements of the operators A and B. 

In the case of the strength matrices, the matrix elements are real for T-even operators 
and image for T-odd operators and thus we easily get 

K p ,K h >0 

K k>l = ~* < 0\[Pk>,X k ]\0 >=4t <ph\Pk>\0><ph\X k ]\Q>, (79) 

ph 

K p ,K h >a 

Vk>l = ~* < Q\[Qk>,Y k ]\0 >= < Ph\Qk>\0 >< ph\Y k ]\0 > (80) 

ph 

where K p and are projections of the momentum of particle and hole states onto quanti- 
zation axis of the system. 

The case of responses is more involved in the sense that matrix elements of the second 
operator in the commutator are transition densities which are generally complex. However, 
the first operator in the commutator still has real (for T-even A) or image (for T-odd A) 
matrix elements and so the averages can be finally reduced to 

K p ,K h >0 

<0|[Q fc ,J a ]|0> = 4i <ph\Qk\0>%<ph\J a \0>, (81) 

ph 

K p ,K h >0 

< 0\[P k ,J a ]\0 > = -4 Y <ph\P k \0>$l<ph\j a \0>, (82) 

ph 

where 3 and 3? result in image and real parts of the transition densities. 
The matrix elements for operator read 

< ph\P k \0 >= i2e ph < ph\Q k \0 > - < ph\Y k \0 > . (83) 
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